home *** CD-ROM | disk | FTP | other *** search
/ Aminet 21 / Aminet 21 (1997)(GTI - Schatztruhe)[!][Oct 1997].iso / Aminet / misc / sci / clustalv.lha / AmigaClustalV / clustalv_sourcecode / trees.c < prev    next >
Encoding:
C/C++ Source or Header  |  1997-07-01  |  16.9 KB  |  679 lines
  1.     /* Phyle of filogenetic tree calculating functions for CLUSTAL V */
  2. #include <stdio.h>
  3. #include <string.h>
  4. #include <stdlib.h>
  5. #include <math.h>
  6. #include "clustalv.h"
  7.  
  8. /*
  9.  *   Prototypes
  10.  */
  11.  
  12. extern void *ckalloc(size_t);
  13. extern int getint(char *, int, int, int);
  14. extern void get_path(char *, char *);
  15. extern FILE * open_output_file(char *, char *, char *, char *);
  16.  
  17.  
  18. void init_trees(void);
  19. void phylogenetic_tree(void);
  20. void bootstrap_tree(void);
  21. void compare_tree(char **, char **, int *, int);
  22. void tree_gap_delete(void); /*flag all positions in alignment that have a gap */
  23. void dna_distance_matrix(FILE *);
  24. void prot_distance_matrix(FILE *);
  25. void nj_tree(char **, FILE *);
  26. void print_tree(char **, FILE *, int *);
  27. void root_tree(char **, int);
  28.  
  29. Boolean transition(int,int);
  30.  
  31. /*
  32.  *   Global variables
  33.  */
  34.  
  35.  
  36. extern double **tmat;     /* general nxn array of reals; allocated from main */
  37.                           /* this is used as a distance matrix */
  38. extern double **smat;
  39. extern Boolean dnaflag;   /* TRUE for DNA seqs; FALSE for proteins */
  40. extern Boolean tossgaps;  /* Ignore places in align. where ANY seq. has a gap*/
  41. extern Boolean kimura;    /* Use correction for multiple substitutions */
  42. extern Boolean empty;     /* any sequences in memory? */
  43. extern Boolean usemenu;   /* interactive (TRUE) or command line (FALSE) */
  44. extern void error(char *,...);  /* error reporting */
  45. extern int nseqs;         /* total no. of seqs. */
  46. extern int *seqlen_array; /* the lengths of the sequences */
  47. extern char **seq_array;   /* the sequences */
  48. extern char **names;       /* the seq. names */
  49. extern char seqname[];        /* name of input file */
  50.  
  51. static double *av;
  52. static int *kill;
  53. static int ran_factor;
  54. int boot_ntrials;            /* number of bootstrap trials */
  55. unsigned int boot_ran_seed;        /* random number generator seed */
  56. static int root_sequence;
  57. static int *boot_positions;
  58. static int *boot_totals;
  59. static char **standard_tree;
  60. static char **sample_tree;
  61. static FILE *phy_tree_file;
  62. static char phy_tree_name[FILENAMELEN+1];
  63. static Boolean verbose;
  64. static char *tree_gaps;    /* array of weights; 1 for use this posn.; 0 don't */
  65.  
  66.  
  67. void init_trees()
  68. {
  69.     register int i;
  70.     
  71.     tree_gaps = (char *)ckalloc( (MAXLEN+1) * sizeof (char) );
  72.         
  73.     boot_positions = (int *)ckalloc( (MAXLEN+1) * sizeof (int) );
  74.     boot_totals    = (int *)ckalloc( (MAXN+1) * sizeof (int) );
  75.  
  76.     kill = (int *) ckalloc( (MAXN+1) * sizeof (int) );
  77.     av   = (double *) ckalloc( (MAXN+1) * sizeof (double)   );
  78.  
  79.     standard_tree = (char **) ckalloc( (MAXN+1) * sizeof (char *) );
  80.     for(i=0; i<MAXN+1; i++) 
  81.         standard_tree[i] = (char *) ckalloc( (MAXN+1) * sizeof(char) );
  82.  
  83.     sample_tree   = (char **) ckalloc( (MAXN+1) * sizeof (char *) );
  84.     for(i=0; i<MAXN+1; i++) 
  85.         sample_tree[i]   = (char *) ckalloc( (MAXN+1) * sizeof(char) );
  86.  
  87.     boot_ntrials  = 1000;
  88.     boot_ran_seed = 111;
  89.         kimura   = FALSE;
  90.         tossgaps = FALSE;    
  91. }
  92.  
  93.  
  94.  
  95. void phylogenetic_tree()
  96. {    char path[FILENAMELEN+1];
  97.     int j;
  98.  
  99.     if(empty) {
  100.         error("You must load an alignment first");
  101.         return;
  102.     }
  103.  
  104.     get_path(seqname,path);
  105.     
  106.     if((phy_tree_file = open_output_file(
  107.         "\nEnter name for tree output file  ",path,
  108.         phy_tree_name,"nj")) == NULL) return;
  109.  
  110.     for(j=1; j<=seqlen_array[1]; ++j) 
  111.         boot_positions[j] = j;        
  112.  
  113.     verbose = TRUE;                  /* Turn on screen output */
  114.     if(dnaflag)
  115.         dna_distance_matrix(phy_tree_file);
  116.     else 
  117.         prot_distance_matrix(phy_tree_file);
  118.  
  119.     verbose = TRUE;              /* Turn on output */
  120.     nj_tree(standard_tree,phy_tree_file);
  121.     fclose(phy_tree_file);    
  122. /*
  123.     print_tree(standard_tree,phy_tree_file);
  124. */
  125.     fprintf(stdout,"\nPhylogenetic tree file created:   [%s]",phy_tree_name);
  126. }
  127.  
  128.  
  129.  
  130.  
  131.  
  132. Boolean transition(int base1, int base2) /* TRUE if transition; else FALSE */
  133. /* 
  134.  
  135.    assumes that the bases of DNA sequences have been translated as
  136.    a,A = 1;   c,C = 2;   g,G = 3;   t,T,u,U = 4;  X or N = 0;  "-" < 0;
  137.  
  138.    A <--> G  and  T <--> C  are transitions;  all others are transversions.
  139.  
  140. */
  141. {
  142.     if( ((base1 == 1) && (base2 == 3)) || ((base1 == 3) && (base2 == 1)) )
  143.         return TRUE;                                     /* A <--> G */
  144.     if( ((base1 == 4) && (base2 == 2)) || ((base1 == 2) && (base2 == 4)) )
  145.         return TRUE;                                     /* T <--> C */
  146.     return FALSE;
  147. }
  148.  
  149.  
  150. void tree_gap_delete()   /* flag all positions in alignment that have a gap */
  151. {              /* in ANY sequence */
  152.     int seqn, posn;
  153.  
  154.     for(posn=1; posn<=seqlen_array[1]; ++posn) {
  155.         tree_gaps[posn] = 0;
  156.          for(seqn=1; seqn<=nseqs; ++seqn)  {
  157.             if(seq_array[seqn][posn] <= 0) {
  158.                tree_gaps[posn] = 1;
  159.                 break;
  160.             }
  161.         }
  162.     }
  163. }
  164.  
  165.  
  166.  
  167. void nj_tree(char **tree_description, FILE *tree)
  168. {
  169.     register int i;
  170.     int l[4],nude,k;
  171.     int nc,mini,minj,j,j1,ii,jj;
  172.     double fnseqs,fnseqs2,sumd;
  173.     double diq,djq,dij,d2r,dr,dio,djo,da;
  174.     double tmin,total,dmin;
  175.     double bi,bj,b1,b2,b3,branch[4];
  176.     int typei,typej,type[4];             /* 0 = node; 1 = OTU */
  177.     
  178.     fnseqs = (double)nseqs;
  179.  
  180. /*********************** First initialisation ***************************/
  181.     
  182.     if(verbose)  {
  183.         fprintf(tree,"\n\n\t\t\tNeighbor-joining Method\n");
  184.         fprintf(tree,"\n Saitou, N. and Nei, M. (1987)");
  185.         fprintf(tree," The Neighbor-joining Method:");
  186.         fprintf(tree,"\n A New Method for Reconstructing Phylogenetic Trees.");
  187.         fprintf(tree,"\n Mol. Biol. Evol., 4(4), 406-425\n");
  188.         fprintf(tree,"\n\n This is an UNROOTED tree\n");
  189.         fprintf(tree,"\n Numbers in parentheses are branch lengths\n\n");
  190.     }    
  191.  
  192.     mini = minj = 0;
  193.  
  194.     for(i=1;i<=nseqs;++i) 
  195.         {
  196.         tmat[i][i] = av[i] = 0.0;
  197.         kill[i] = 0;
  198.         }
  199.  
  200. /*********************** Enter The Main Cycle ***************************/
  201.  
  202.  /*    for(nc=1; nc<=(nseqs-3); ++nc) {  */                /**start main cycle**/
  203.     for(nc=1; nc<=(nseqs-3); ++nc) {
  204.         sumd = 0.0;
  205.         for(j=2; j<=nseqs; ++j)
  206.             for(i=1; i<j; ++i) {
  207.                 tmat[j][i] = tmat[i][j];
  208.                 sumd = sumd + tmat[i][j];
  209.                 smat[i][j] = smat[j][i] = 0.0;
  210.             }
  211.  
  212.         tmin = 99999.0;
  213.  
  214. /*.................compute SMATij values and find the smallest one ........*/
  215.  
  216.         for(jj=2; jj<=nseqs; ++jj) 
  217.             if(kill[jj] != 1) 
  218.                 for(ii=1; ii<jj; ++ii)
  219.                     if(kill[ii] != 1) {
  220.                         diq = djq = 0.0;
  221.  
  222.                         for(i=1; i<=nseqs; ++i) {
  223.                             diq = diq + tmat[i][ii];
  224.                             djq = djq + tmat[i][jj];
  225.                         }
  226.  
  227.                         dij = tmat[ii][jj];
  228.                         d2r = diq + djq - (2.0*dij);
  229.                         dr  = sumd - dij -d2r;
  230.                         fnseqs2 = fnseqs - 2.0;
  231.                             total= d2r+ fnseqs2*dij +dr*2.0;
  232.                         total= total / (2.0*fnseqs2);
  233.                         smat[ii][jj] = total;
  234.  
  235.                         if(total < tmin) {
  236.                             tmin = total;
  237.                             mini = ii;
  238.                             minj = jj;
  239.                         }
  240.                     }
  241.         
  242.  
  243. /*.................compute branch lengths and print the results ........*/
  244.  
  245.  
  246.         dio = djo = 0.0;
  247.         for(i=1; i<=nseqs; ++i) {
  248.             dio = dio + tmat[i][mini];
  249.             djo = djo + tmat[i][minj];
  250.         }
  251.  
  252.         dmin = tmat[mini][minj];
  253.         dio = (dio - dmin) / fnseqs2;
  254.         djo = (djo - dmin) / fnseqs2;
  255.         bi = (dmin + dio - djo) * 0.5;
  256.         bj = dmin - bi;
  257.         bi = bi - av[mini];
  258.         bj = bj - av[minj];
  259.  
  260.         if( av[mini] > 0.0 )
  261.             typei = 0;
  262.         else
  263.             typei = 1;
  264.         if( av[minj] > 0.0 )
  265.             typej = 0;
  266.         else
  267.             typej = 1;
  268.  
  269.         if(verbose) {
  270.              fprintf(tree,"\n Cycle%4d     = ",nc);
  271.             if(typei == 0)
  272.             fprintf(tree,"Node:%4d (%9.5f) joins ",mini,bi);
  273.             else
  274.              fprintf(tree," SEQ:%4d (%9.5f) joins ",mini,bi);
  275.             if(typej == 0) 
  276.             fprintf(tree,"Node:%4d (%9.5f)",minj,bj);
  277.             else
  278.             fprintf(tree," SEQ:%4d (%9.5f)",minj,bj);
  279.             fprintf(tree,"\n");
  280.         }
  281.  
  282.         for(i=1; i<=nseqs; i++)
  283.             tree_description[nc][i] = 0;
  284.  
  285.              if(typei == 0) { 
  286.             for(i=nc-1; i>=1; i--)
  287.                 if(tree_description[i][mini] == 1) {
  288.                     for(j=1; j<=nseqs; j++)  
  289.                          if(tree_description[i][j] == 1)
  290.                             tree_description[nc][j] = 1;
  291.                     break;
  292.                 }
  293.         }
  294.         else
  295.             tree_description[nc][mini] = 1;
  296.  
  297.         if(typej == 0) {
  298.             for(i=nc-1; i>=1; i--) 
  299.                 if(tree_description[i][minj] == 1) {
  300.                     for(j=1; j<=nseqs; j++)  
  301.                          if(tree_description[i][j] == 1)
  302.                             tree_description[nc][j] = 1;
  303.                     break;
  304.                 }
  305.         }
  306.         else
  307.             tree_description[nc][minj] = 1;
  308.             
  309.         if(dmin <= 0.0) dmin = 0.0001;
  310.         av[mini] = dmin * 0.5;
  311.  
  312. /*........................Re-initialisation................................*/
  313.  
  314.         fnseqs = fnseqs - 1.0;
  315.         kill[minj] = 1;
  316.  
  317.         for(j=1; j<=nseqs; ++j) 
  318.             if( kill[j] != 1 ) {
  319.                 da = ( tmat[mini][j] + tmat[minj][j] ) * 0.5;
  320.                 if( (mini - j) < 0 ) 
  321.                     tmat[mini][j] = da;
  322.                 if( (mini - j) > 0)
  323.                     tmat[j][mini] = da;
  324.             }
  325.  
  326.         for(j=1; j<=nseqs; ++j)
  327.             tmat[minj][j] = tmat[j][minj] = 0.0;
  328.  
  329.  
  330. /****/    }                        /**end main cycle**/
  331.  
  332. /******************************Last Cycle (3 Seqs. left)********************/
  333.  
  334.     nude = 1;
  335.  
  336.     for(i=1; i<=nseqs; ++i)
  337.         if( kill[i] != 1 ) {
  338.             l[nude] = i;
  339.             nude = nude + 1;
  340.         }
  341.  
  342.     b1 = (tmat[l[1]][l[2]] + tmat[l[1]][l[3]] - tmat[l[2]][l[3]]) * 0.5;
  343.     b2 =  tmat[l[1]][l[2]] - b1;
  344.     b3 =  tmat[l[1]][l[3]] - b1;
  345.     branch[1] = b1 - av[l[1]];
  346.     branch[2] = b2 - av[l[2]];
  347.     branch[3] = b3 - av[l[3]];
  348.  
  349.  
  350.     for(i=1; i<=nseqs; i++)
  351.         tree_description[nseqs-2][i] = 0;
  352.  
  353.     if(verbose)
  354.         fprintf(tree,"\n Cycle%4d (Last cycle, trichotomy):\n",nc);
  355.  
  356.     for(i=1; i<=3; ++i) {
  357.        if( av[l[i]] > 0.0) {
  358.               if(verbose)
  359.                   fprintf(tree,"\n\t\t Node:%4d (%9.5f) ",l[i],branch[i]);
  360.         for(k=nseqs-3; k>=1; k--)
  361.             if(tree_description[k][l[i]] == 1) {
  362.                 for(j=1; j<=nseqs; j++)
  363.                      if(tree_description[k][j] == 1)
  364.                         tree_description[nseqs-2][j] = i;
  365.                 break;
  366.             }
  367.        }
  368.        else  {
  369.               if(verbose)
  370.                fprintf(tree,"\n\t\t  SEQ:%4d (%9.5f) ",l[i],branch[i]);
  371.         tree_description[nseqs-2][l[i]] = i;
  372.        }
  373.        if(i < 3) {
  374.               if(verbose)
  375.                 fprintf(tree,"joins");
  376.        }
  377.     }
  378.  
  379.     if(verbose)
  380.         fprintf(tree,"\n");
  381.  
  382. }
  383.  
  384.  
  385.  
  386.  
  387. void bootstrap_tree()
  388. {
  389.     int i,j,ranno;
  390.     int k,l,m,p;
  391.     unsigned int num;
  392.     char lin2[MAXLINE],path[MAXLINE+1];
  393.  
  394.     if(empty) {
  395.         error("You must load an alignment first");
  396.         return;
  397.     }
  398.  
  399.     get_path(seqname, path);
  400.     
  401.     if((phy_tree_file = open_output_file(
  402.         "\nEnter name for bootstrap output file  ",path,
  403.         phy_tree_name,"njb")) == NULL) return;
  404.  
  405.     for(i=0;i<MAXN+1;i++)
  406.         boot_totals[i]=0;
  407.         
  408.     for(j=1; j<=seqlen_array[1]; ++j)  /* First select all positions for */
  409.         boot_positions[j] = j;       /* the "standard" tree */
  410.  
  411.     verbose = TRUE;                   /* Turn on screen output */
  412.     if(dnaflag)
  413.         dna_distance_matrix(phy_tree_file);
  414.     else 
  415.         prot_distance_matrix(phy_tree_file);
  416.  
  417.     verbose = TRUE;                  /* Turn on screen output */
  418.     nj_tree(standard_tree, phy_tree_file);    /* compute the standard tree */
  419.  
  420.     fprintf(phy_tree_file,"\n\n\t\t\tBootstrap Confidence Limits\n\n");
  421.  
  422.     ran_factor = RAND_MAX / seqlen_array[1];
  423.  
  424.     if(usemenu) 
  425.            boot_ran_seed = 
  426. getint("\n\nEnter seed no. for random number generator ",1,1000,boot_ran_seed);
  427.  
  428.     srand(boot_ran_seed);
  429.     fprintf(phy_tree_file,"\n Random number generator seed = %7u\n",
  430.     boot_ran_seed);
  431.  
  432.     if(usemenu) 
  433.           boot_ntrials = 
  434. getint("\n\nEnter number of bootstrap trials ",1,10000,boot_ntrials);
  435.  
  436.       fprintf(phy_tree_file,"\n Number of bootstrap trials   = %7d\n",
  437.     boot_ntrials);
  438.  
  439.     fprintf(phy_tree_file,
  440.     "\n\n Diagrammatic representation of the above tree: \n");
  441.     fprintf(phy_tree_file,"\n Each row represents 1 tree cycle;");
  442.     fprintf(phy_tree_file," defining 2 groups.\n");
  443.     fprintf(phy_tree_file,"\n Each column is 1 sequence; ");
  444.     fprintf(phy_tree_file,"the stars in each line show 1 group; ");
  445.     fprintf(phy_tree_file,"\n the dots show the other\n");
  446.     fprintf(phy_tree_file,"\n Numbers show occurences in bootstrap samples.");
  447. /*
  448.     print_tree(standard_tree, phy_tree_file, boot_totals);
  449. */
  450.     verbose = FALSE;                   /* Turn OFF screen output */
  451.  
  452.     fprintf(stdout,"\n\nEach dot represents 10 trials\n\n");
  453.     for(i=1; i<=boot_ntrials; ++i) {
  454.         for(j=1; j<=seqlen_array[1]; ++j) {      /* select alignment */
  455.             ranno = ( rand() / ran_factor ) + 1; /* positions for */
  456.             boot_positions[j] = ranno;       /* bootstrap sample */
  457.         }
  458.         if(dnaflag)
  459.             dna_distance_matrix(phy_tree_file);
  460.         else
  461.             prot_distance_matrix(phy_tree_file);
  462.         nj_tree(sample_tree, phy_tree_file); /* compute 1 sample tree */
  463.         compare_tree(standard_tree, sample_tree, boot_totals, nseqs);
  464.         if(i % 10  == 0) fprintf(stdout,".");
  465.         if(i % 100 == 0) fprintf(stdout,"\n");
  466.     }
  467.  
  468. /*
  469.     fprintf(phy_tree_file,"\n\n Bootstrap totals for each group\n");
  470. */
  471.     print_tree(standard_tree, phy_tree_file, boot_totals);
  472.  
  473.     fclose(phy_tree_file);
  474.  
  475.     fprintf(stdout,"\n\nBootstrap output file completed       [%s]"
  476.         ,phy_tree_name);
  477. }
  478.  
  479.  
  480. void compare_tree(char **tree1, char **tree2, int *hits, int n)
  481. {    
  482.     int i,j,k,l;
  483.     int nhits1, nhits2;
  484.  
  485.     for(i=1; i<=n-3; i++)  {
  486.         for(j=1; j<=n-3; j++)  {
  487.             nhits1 = 0;
  488.             nhits2 = 0;
  489.             for(k=1; k<=n; k++) {
  490.                 if(tree1[i][k] == tree2[j][k]) nhits1++;
  491.                 if(tree1[i][k] != tree2[j][k]) nhits2++;
  492.             }
  493.             if((nhits1 == nseqs) || (nhits2 == nseqs)) hits[i]++;
  494.         }
  495.     }
  496. }
  497.  
  498.  
  499.  
  500.  
  501. void print_tree(char **tree_description, FILE *tree, int *totals)
  502. {
  503.     int row,col;
  504.  
  505.     fprintf(tree,"\n");
  506.  
  507.     for(row=1; row<=nseqs-3; row++)  {
  508.         fprintf(tree," \n");
  509.         for(col=1; col<=nseqs; col++) { 
  510.             if(tree_description[row][col] == 0)
  511.                 fprintf(tree,"*");
  512.             else
  513.                 fprintf(tree,".");
  514.         }
  515.         if(totals[row] > 0)
  516.             fprintf(tree,"%7d",totals[row]);
  517.     }
  518.     fprintf(tree," \n");
  519.     for(col=1; col<=nseqs; col++) 
  520.         fprintf(tree,"%1d",tree_description[nseqs-2][col]);
  521.     fprintf(tree,"\n");
  522. }
  523.  
  524.  
  525.  
  526. void dna_distance_matrix(FILE *tree)
  527. {   
  528.     int m,n,j,i;
  529.     int res1, res2;
  530.     double p,q,e,a,b,k;    
  531.  
  532.     tree_gap_delete();  /* flag positions with gaps (tree_gaps[i] = 1 ) */
  533.     
  534.     if(verbose) {
  535.         fprintf(tree,"\n");
  536.         fprintf(tree,"\n DIST   = percentage divergence (/100)");
  537.         fprintf(tree,"\n p      = rate of transition (A <-> G; C <-> T)");
  538.         fprintf(tree,"\n q      = rate of transversion");
  539.         fprintf(tree,"\n Length = number of sites used in comparison");
  540.         fprintf(tree,"\n");
  541.         if(tossgaps) {
  542.         fprintf(tree,"\n All sites with gaps (in any sequence) deleted!");
  543.         fprintf(tree,"\n");
  544.         }
  545.         if(kimura) {
  546.         fprintf(tree,"\n Distances corrected by Kimura's 2 parameter model:");
  547.         fprintf(tree,"\n\n Kimura, M. (1980)");
  548.         fprintf(tree," A simple method for estimating evolutionary ");
  549.         fprintf(tree,"rates of base");
  550.         fprintf(tree,"\n substitutions through comparative studies of ");
  551.         fprintf(tree,"nucleotide sequences.");
  552.         fprintf(tree,"\n J. Mol. Evol., 16, 111-120.");
  553.         fprintf(tree,"\n\n");
  554.         }
  555.     }
  556.  
  557.     for(m=1;   m<nseqs;  ++m)     /* for every pair of sequence */
  558.     for(n=m+1; n<=nseqs; ++n) {
  559.         p = q = e = 0.0;
  560.         tmat[m][n] = tmat[n][m] = 0.0;
  561.         for(i=1; i<=seqlen_array[1]; ++i) {
  562.             j = boot_positions[i];
  563.                         if(tossgaps && (tree_gaps[j] > 0) ) 
  564.                 goto skip;          /* gap position */
  565.             res1 = seq_array[m][j];
  566.             res2 = seq_array[n][j];
  567.             if( (res1 < 1) || (res2 < 1) )  
  568.                 goto skip;          /* gap in a seq*/
  569.             e = e + 1.0;
  570.                         if(res1 != res2) {
  571.                 if(transition(res1,res2))
  572.                     p = p + 1.0;
  573.                 else
  574.                     q = q + 1.0;
  575.             }
  576.                 skip:;
  577.         }
  578.  
  579.  
  580.     /* Kimura's 2 parameter correction for multiple substitutions */
  581.  
  582.         if(!kimura) {
  583.             k = (p+q)/e;
  584.             if(p > 0.0)
  585.                 p = p/e;
  586.             else
  587.                 p = 0.0;
  588.             if(q > 0.0)
  589.                 q = q/e;
  590.             else
  591.                 q = 0.0;
  592.             tmat[m][n] = tmat[n][m] = k;
  593.             if(verbose)                    /* if screen output */
  594.                 fprintf(tree,        
  595.           "%4d vs.%4d:  DIST = %7.4f; p = %6.4f; q = %6.4f; length = %6.0f\n"
  596.                              ,m,n,k,p,q,e);
  597.         }
  598.         else {
  599.             if(p > 0.0)
  600.                 p = p/e;
  601.             else
  602.                 p = 0.0;
  603.             if(q > 0.0)
  604.                 q = q/e;
  605.             else
  606.                 q = 0.0;
  607.             a = 1.0/(1.0-2.0*p-q);
  608.             b = 1.0/(1.0-2.0*q);
  609.             k = 0.5*log(a) + 0.25*log(b);
  610.             tmat[m][n] = tmat[n][m] = k;
  611.             if(verbose)                      /* if screen output */
  612.                    fprintf(tree,
  613.              "%4d vs.%4d:  DIST = %7.4f; p = %6.4f; q = %6.4f; length = %6.0f\n"
  614.                             ,m,n,k,p,q,e);
  615.  
  616.         }
  617.     }
  618. }
  619.  
  620.  
  621.  
  622. void prot_distance_matrix(FILE *tree)
  623. {
  624.     int m,n,j,i;
  625.     int res1, res2;
  626.     double p,e,a,b,k;    
  627.  
  628.     tree_gap_delete();  /* flag positions with gaps (tree_gaps[i] = 1 ) */
  629.     
  630.     if(verbose) {
  631.         fprintf(tree,"\n");
  632.         fprintf(tree,"\n DIST   = percentage divergence (/100)");
  633.         fprintf(tree,"\n Length = number of sites used in comparison");
  634.         fprintf(tree,"\n\n");
  635.             if(tossgaps) {
  636.             fprintf(tree,"\n All sites with gaps (in any sequence) deleted");
  637.             fprintf(tree,"\n");
  638.         }
  639.             if(kimura) {
  640.             fprintf(tree,"\n Distances corrected by Kimura's empirical method:");
  641.             fprintf(tree,"\n\n Kimura, M. (1983)");
  642.             fprintf(tree," The Neutral Theory of Molecular Evolution.");
  643.             fprintf(tree,"\n Cambridge University Press, Cambridge, England.");
  644.             fprintf(tree,"\n\n");
  645.             }
  646.     }
  647.  
  648.     for(m=1;   m<nseqs;  ++m)     /* for every pair of sequence */
  649.     for(n=m+1; n<=nseqs; ++n) {
  650.         p = e = 0.0;
  651.         tmat[m][n] = tmat[n][m] = 0.0;
  652.         for(i=1; i<=seqlen_array[1]; ++i) {
  653.             j = boot_positions[i];
  654.                     if(tossgaps && (tree_gaps[j] > 0) ) goto skip; /* gap position */
  655.             res1 = seq_array[m][j];
  656.             res2 = seq_array[n][j];
  657.             if( (res1 < 1) || (res2 < 1) )  goto skip;   /* gap in a seq*/
  658.             e = e + 1.0;
  659.                         if(res1 != res2) p = p + 1.0;
  660.                 skip:;
  661.         }
  662.  
  663.         if(p <= 0.0) 
  664.             k = 0.0;
  665.         else
  666.             k = p/e;
  667.  
  668.         if(kimura) 
  669.             if(k > 0.0) k = - log(1.0 - k - (k * k/5.0) );
  670.  
  671.         tmat[m][n] = tmat[n][m] = k;
  672.             if(verbose)                    /* if screen output */
  673.             fprintf(tree,        
  674.                       "%4d vs.%4d  DIST = %6.4f;  length = %6.0f\n",m,n,k,e);
  675.     }
  676. }
  677.  
  678.  
  679.